Collisional statistics of the hard-sphere gas 



Paolo Visco, 12 Frederic van Wijland, 3 and Emmanuel Trizac 1 

' Universite Paris-Sud, LPTMS, UMR 8626, Orsay Cedex, F -91405 and CNRS, Orsay, F-91405 
2 SUPA, School of Physics, University of Edinburgh, Mayfield Road, Edinburgh, EH9 3JZ, UK 
3 Laboratoire Matiere et Systemes Complexes (CNRS UMR 7057), Universite Denis Diderot (Paris VII), 
10 rue Alice Domon et Leonie Duquet, 75205 Paris cedex 13, France 
(Dated: March 9, 2008) 

We investigate the probability distribution function of the free flight time and of the number of collisions in a 
hard sphere gas at equilibrium. At variance with naive expectation, the latter quantity does not follow Poissonian 
statistics, even in the dilute limit which is the focus of the present analysis. The corresponding deviations are 
addressed both numerically and analytically. In writing an equation for the generating function of the cumulants 
of the number of collisions, we came across a perfect mapping between our problem and a previously introduced 
model: the probabilistic ballistic annihilation process [Coppex et al, Phys. Rev. E 69 11303 (2004)]. We 
exploit this analogy to construct a Monte-Carlo like algorithm able to investigate the asymptotically large time 
behavior of the collisional statistics within a reasonable computational time. In addition, our predictions are 
confronted against the results of Molecular Dynamics simulations and Direct Simulation Monte Carlo technique. 
An excellent agreement is reported. 

PACS numbers: 05.40.-a, 51.10.+y, 02.50.-r 



I. INTRODUCTION 

Despite the great interest that the hard sphere gas has trig- 
gered since the early days of statistical physics, there are still, 
to this day, simple questions which have not been deeply in- 
vestigated yet. In this paper we address one such question: 
what is the probability that a tagged particle suffers a given 
number of collisions in a time tl Apart from the simplicity of 
the question, relevant in its own right, the collisional statistics 
of the hard sphere gas turns out to be useful in several pur- 
poses, such as the estimation of transport coefficients or the 
characterization of the dissipation in the closely related gran- 
ular gases. 

In a low density hard sphere gas at equilibrium, the veloci- 
ties of two colliding particles are uncorrelated just before the 
collision: this is the molecular chaos (strosszahlansatz) state- 
ment. This remark could naively lead to the conclusion that 
collisions are uncorrelated random events, implying therefore 
that the number of collisions is simply a Poisson random vari- 
able. However, even if molecular chaos is exactly verified, the 
collision number is not a Poisson random variable. The non- 
Poissonian nature of collisional statistics has already been no- 
ticed in the literature (see e.g. [1]) but to our knowledge, has 
resisted analytical investigations. It is our purpose here to fill 
this gap. 

The reason for the non-Poissonian behavior alluded to 
above is that the probability for a collision to take place de- 
pends on the scattering cross section of the colliding pair 
( which depends on the relative velocity g of the pair). Namely, 
for the hard sphere gas, the probability of a collision behaves 
roughly as g • <rO(g ■ er), where <r is a unitary vector joining 
the centers of mass of the two particles at contact (directed 
along the apse line), and O is the Heaviside step function. 
For a gas made up of particles interacting through a poten- 
tial other than the hard-sphere potential, the probability of 
having a "collision" is of course different, which affects the 
distribution of the number of collisions. Of particular inter- 



est is the gas of particle interacting through a pair potential 
V(r) ~ l/r 2 ^ -1 ), where d is the space dimension. Such 
particles are known as Maxwell molecules yfl and lead to a 
velocity independent probability of having a collision. Within 
this model, it then appears that the collisions are truly uncor- 
related random events, and hence that the number of collision 
is a Poisson random variable. 



In the next section, the free flight time distribution of a hard 
sphere gas is addressed, together with the distribution of free 
path lengths. These results pertain to AT — collision proper- 
ties. In sectionUnl the analysis is extended to consider the full 
probability P(Af, t) encoding the number of collisions statis- 
tics. The large time behavior and a complementary perturba- 
tive treatment are worked out analytically. Explicit expres- 
sions are obtained for the cumulants of the number of colli- 
sions J\f. In order to put the theoretical predictions to the test, 
three types of numerical simulations are performed. The first 
two, Molecular Dynamics and Direct Simulation Monte Carlo 
(DSMC) are routinely employed in the field and beyond. The 
third type, of the Monte Carlo family, is discussed in llVl and 
specifically devised to solve the particular problem at hand. It 
relies on a reinterpretation of the eigenvalue equation for the 
generating function for the cumulants of A/, in terms of a pop- 
ulation dynamics with cloning and annihilation. A Markov 
chain with the desired properties in then constructed, which 
allows for a direct measure of several key quantities involved 
in the analytical treatment. In doing so, we uncover a fruit- 
ful mapping with the probabilistic ballistic annihilation model 
proposed in Ref [3]. The three numerical methods provide re- 
sults that are in excellent agreement with the analytical predic- 
tions. Preliminary accounts of part of this work has appeared 
in ||4j], where the numerical aspect of the work was limited to 
its Molecular Dynamics content, and where the theory is re- 
stricted to the zeroth order of the treatment put forward here. 



2 



H. FREE FLIGHTS TIME DISTRIBUTION 
A. General results 

We consider a hard sphere gas in d dimensions composed 
of N particles of equal mass m and equal diameter a. The 
gas is thermalized at some temperature To in a homogeneous 
state of (constant) density p, so that the one point distribution 
function of the system is a Gaussian: 



<Kr,v) 



(v)=p(27rT )- d / 2 exp 
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2T, 



(1) 



If one decides to follow the evolution of one particle among 
the N, then the velocity probability distribution function (pdf) 
of that particle will evolve according to the homogeneous lin- 
ear Boltzmann equation: 
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where I — (<r d ~ 1 px) 1 is a length proportional to the mean 
free path denoted / below, and V12 = vi — V2 . The fac- 
tor x is the pair correlation function at contact and is the so 
called Enskog correction factor J5|]. In the dilute limit where 
the molecular chaos assumption is justified, it tends to unity. 
Moreover, the primed integral means that the integration has 
to be performed in the domain for which V12 • o > 0. Finally, 
we are employing the short-hand notation /** = /(v**,t), 
where the two star superscript refers to the precollisional ve- 
locity: 



Vl - (V12 ■ <r)(T , 



v 2 + (v 12 ■ <t)ct . (3) 



In this framework the evolution of our tagged particle is ex- 
actly a Markov process. The probability of hopping from a 
state with velocity v to another state of velocity in a narrow 
interval dv' around v' and in a time interval di is given by 
W(v'\v)dv'dt, where 



1 



W(V\w) = ^|v'-v| 2 - d (27rT )- 1 / 2 exp 
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is the transition rate density per unit time (see e.g. Ig]). For 
any general Markov process, the probability of leaving a given 
configuration T in a time interval dt is r(T)dt, where 



r(r) = J dT'W{T'\T) 



(5) 



This is simply the loss term of the Linear Boltzmann equation 
(O, and it reads for the hard sphere gas: 
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where \Fx is the confluent hyper-geometric function of the 
first kind Q7|] , and u> is the collision frequency of the gas: 



u) = I dv0(v)r(w) = 
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(7) 



From the Markovian property, it follows that the probability 
of having a given velocity v for a given time t and making a 
collision exactly at time t is exponential: 

Pfft(*|v) = r(v) exp (-r(v)t) . (8) 

The probability density of the free flight time follows: 



FFT 



(t) 
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Here ^ co /i(v) is the velocity pdf of the colliding particles (i.e. 
obtained sampling the velocity of the particle only on colli- 
sion): 



r(v) 



(10) 



where the normalization factor uj is exactly the collision fre- 
quency of the gas. The free flight time distribution for hard 
spheres was already investigated in [8], but the authors used 
the equilibrium velocity distribution (f> instead of the on col- 
lision velocity distribution 4> co u as a weight in (0, and hence 
obtained an incorrect result (this issue is further commented 
in Appendix [A}. Two measurements of the free flight time 
distribution in an event driven molecular dynamic simulation 
of hard disks (see section [IV] for more details) and in a Direct 
Simulation Monte Carlo (DSMC) IJ, are shown in Fig. [Q 
together with the expression of Eq. The agreement is ex- 
cellent. The molecular dynamics simulations are performed 
at a low density (pa 2 = 10~ 2 ) to ensure that molecular chaos 
holds, whereas the DSMC technique relies by construction on 
the molecular chaos approximation, and can therefore be con- 
sidered as providing a p — > benchmark. Figure Q] also pro- 
vides a comparison with the exponential law that would hold 
for a constant r(v), as is the case for Maxwell molecules. The 
latter expectation [oj exp(— cut)] is seen to hold at early times, 
but significantly fails at large t, a time regime that we investi- 
gate in more details in the next section. 



B. Large time and high dimension analysis 

Apart from simplified models (see AppendixlBl. it does not 
seem possible to obtain a simple exact analytical expression 
for the free flight time distribution, mainly because of the par- 
ticular behavior of the hyper-geometric function involved in 
the expression of the collision rate r(v). In order to get a sim- 
pler approximated result, one has to look at limiting cases. 
For large times, the integral involved in Eq. <j9j can be esti- 
mated through the saddle-point approximation (a.k.a. method 
of steepest descent, see e.g. 111011 ). yielding: 
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FIG. 1: (Color online) Free flights time distribution of a hard disc 
gas. The circles and squares correspond to the results of Molec- 
ular Dynamics (MD) and Direct Simulation Monte-Carlo (DSMC) 
methods. MD simulations in all this work are shown at density 
p = 0.01a -2 . The full line is the numerical integration of Eq. ((9}. 
The dot-dashed line shows the asymptotic large time behavior of Eq. 
{TTJ, and the dashed line is the exponential distribution, associated to 
the Poisson distribution. Time is measured in units of the mean free 
time {uj — 1). 



The above function is plotted in Fig. Q] and seen to be in 
very good agreement with our numerical data provided t is 
large enough. Note that expression ( fTTT l does not provide a 
normalized probability (a feature visible on Fig. [T), but only 
an asymptotic expansion. 

Besides, a way to get an explicit expression for the collision 
rate r(v) is to investigate the infinite dimension limit. In this 
case, the saddle-point approximation yields: 



r \ UJ I 
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(12) 



A detailed derivation of this result is given in AppendixICl 



C. Free path length distribution 

From the knowledge of the free flight time distribution, it 
is straightforward to compute the distribution of the free path 
length. In fact, for a particle with a velocity v, the length cov- 
ered in a time t is x = \v\t = vt. Therefore, the conditional 
distribution of free path length is obtained as: 
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The free path length distribution can be obtained averaging 
the above expression with the weight 4> co u: 



Pfpl(x) = J Av(j) co u{-v)P F p L {x\v) 



(14) 
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FIG. 2: (Color online) Plot of the free path length distribution from 
MD simulations (at density p — 0.01<t -2 ) and from a numerical 
integration of Eq. d 1 4b . In the inset, the dashed line is the dominant 
exponential term, Eq. dl6t . 



Although the conditional probability of free path length is 
very similar to the conditional probability of free flight time, 
the step of averaging over the collisional velocity distribution 
may drastically change the shape of the distribution. A case 
worthy of attention is discussed in Appendix iDl for Maxwell 
molecules: the free flight time distribution is a pure expo- 
nential, while the free path length distribution has a stretched 
exponential decay. In the case of hard particles, the domi- 
nant large path behavior remains exponential, as for the free 
flight time. Nevertheless, the leading exponential prefactor 
changes. The free flight time distribution is dominated by 
e -r(o)t^ wnere r (rj) is the minimum of the function r(v). In 
the case of free path length, the dominant exponential behav- 
ior is s Ax , where A is the minimum of r(v)/v, which is ob- 
tained in the limit v — > oo: 



A = lim 



r(v) 



1 

771 



7T 2 



(15) 



where I = u>/(v) is the mean free path. This yields: 



FPL 



(x) 



e V21 



(16) 



(here the ~ should be understood as an equivalence between 
the logarithms in the large x limit). Unfortunately, this expres- 
sion only gives the behavior of the dominating exponential 
term, which becomes visible at a length scales significantly 
larger than the mean free path. Hence, if one analyzes the 
result of MD simulations as in Fig. |2j one sees that there 
are sub-leading (algebraic) terms in the asymptotic behavior, 
which still play a role and that are responsible for the mis- 
match between Eq. dT6b and MD data in Fig. [2] A similar 
feature holds for the free flight time distribution: upon ne- 
glecting the algebraic sub-dominant correction in Eq. ( fTTI ). 
the agreement displayed at late times in the inset of Fig. Q] 
would be spoiled. 
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III. THE NUMBER OF COLLISIONS 



B. Large time behavior 



A. Uncorrected approach 

After having shown how the probability of the time between 
two subsequent collisions behaves, our interest goes to the 
probability of the sum of many such times, and in particu- 
lar, the probability of having Af collisions in a time t. The 
time interval t in which the J\f collisions take place can be 
decomposed in a sum of Af + 1 (correlated) times: 



t = ti n + ti + • ■ ■ + tj^-l + tf . 



(17) 



We are considering that a given particle begins a trajectory 
with a given initial velocity, and then waits a time t in before 
it makes a collision and gets a new velocity. During a time 
ti, it then flies freely, before colliding again etc. The succes- 
sive behavior of the particle's trajectory will be a sequence of 
collisions, spaced by other free flight times i 2 , t%, . . . , tjs/-\, 
until the last 7V-th collision, which will be followed by a final 
time interval t f which is not ended with a collision. The free 
flight times t\ to are of course distributed following the 
probability Pfft given by Eq. We will denote by Pi n 
and Pf the pdf of t{ n and tf, although their precise statistics 
are irrelevant for our purposes, as becomes clear below. 

If one would assume these Af + 1 times to be uncorrected, 
then the probability P(AT, t) of having Af collisions in a time 
t could be deduced from the knowledge of the free flight time 
probability as follows. We first introduce the Laplace trans- 
form: 



P(z) = / At e- zt P(t) 



(18) 



where P can be either Pfft, Pin, or Pf. The transform of 
P(Af, t) can hence be expressed as: 



P(Af,z) = P rn {z) Pfft^Y- 1 P f {z) . 



(19) 



Of course, since Pfft cannot be computed analytically, a 
closed form expression for P(Af, z) is not available. One 
interesting limit to analyze is that of large times. To get an 
asymptotic form of the following integral 

P(Af,t) = Jdz t zt P m {z) PpFTiz)^- 1 P f (z) , (20) 

we make the remark that the number of collisions increases 
with time linearly on average, and we define a time intensive 
counterpart to Af, a fluctuating collision rate, which we denote 
by n = Af/t. The integral of Eq. ( |20b then reads: 



P(Af,t)= / dz exp(t(z + nlii(P P FT(z)))+0(l) 



(21) 

Finally, when t — > oo, the above probability behaves as: 

P(Af,t) ~e ir W* , (22) 

where 7r(n) is a large deviation function, and is related to 
Pp ft through the following relation: 

7r(n) = min(z + n \ti(Pfft(z))) • (23) 



The result of the previous subsection leads to an approx- 
imate evaluation of the large deviation function of the num- 
ber of collisions that neglects temporal correlations. We will 
now investigate this large deviation function keeping into ac- 
count these correlations. Let us define the joint probability 
f(v,Af, t) of having a velocity v and having suffered Af col- 
lisions up until time t. In a homogeneous state, the time evolu- 
tion of the above defined probability is governed by a slightly 
different form of the linear Boltzmann equation: 



d t f(v x ,Af,t) 



dv? 



d<r(vi2 ■ <x) x 



[/«,-V-l,t)^(v5*)-/(vi,jV,t)0(v 2 )] . (24) 

The function </>(v) still represents the one point velocity pdf of 
the gas, which is in equilibrium. Moreover, since the particle 
whose collisions we are counting has the same velocity pdf, 
we enforce that: 



£/(v,A0 = *(v). 
It is useful to introduce the generating function of / as: 



oo 

^ — AW . 



f(v,Af,t). 
It can be seen that / evolves according to 



/(v,A,t)= 



(25) 



(26) 



dv 2 / d<r6(vi2 • <x)(vi 2 • <x)x 



fit/(vi,A,t) 

[e- A /(vr,A,O0(vr)-/(vi,A,t)0(v 2 )] . (27) 

The large time behavior of the solution of the above equation 
is dominated by the largest eigenvalue /i(A) of the evolution 
operator: 

M (A)/(v 1 ,A) = L A /(v 1 ,A) = 

J dv 2 J der9(vi 2 • <r)(vi 2 • <x)x 

[e- A /(vr,A)0(vr)-/(vi,A)0(v 2 )] , (28) 

where / denotes the eigenfunction of L\ associated with [i\ 

f(v,\t) oc e^ A ) { /(v,A). (29) 

Moreover, since P(A) ~ e^)*, one sees that /u(A) is propor- 
tional to the cumulant generating function: 



(Af p )c = t(-iy 



(30) 



A=0 



Furthermore, P{Af) ~ J dAe AA ^P(A), and hence, applying 
the saddle point method for t — > oo, one has that 

P(Af,t) ~ exp(tn(n)) , (31) 

where n = Af/t and the large deviation function tt is related 
to /z through a Legendre transform: 

7r(n) = min(/i(A) + An) . (32) 

A 
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1. Large X behavior 



Before trying to solve Eq. d28l > extending the methods of 
kinetic theory, we shall first provide exact results, which can 
be extracted from the analysis of asymptotically large values 
of A. The evolution operator appearing in Eq. ( 1281 can be cast 
in the form: 



The result obtained in Eq. (1391 1 is compatible with the large 
time behavior embodied in Eq. ( fTTT i: the time derivative of the 
probability of having Af = collisions is of course (minus) 
the free flight time distribution function. We emphasize that 
both results ( |39l and (fTTT l hold for large times. In addition, 
we note that the only space dimension dependence involved is 
through the collision frequency ui [see Eq. (0]. 



L\ — Ln 



(33) 



where 



£o/(vi,A) 



dv 2 / d<r(vi2 • o-)/(vi, A)</>(v 2 ) = 
-r(«i)/(vi,A) , (34a) 



L 1 /(v 1 ,A) = J dv 2 / d*(vi2-*)/(vr,A)0(vr) • 

(34b) 

For large values of A, the coefficient e _A plays the role of a 
small parameter, and therefore the eigenvalue equation (|28l 
can be solved in perturbation theory. We will therefore try to 
get an expression of /i(A) for large A of the form: 



M (A) = /z^ + e" A /i (1) + <^(e~ 2A ) . 



(35) 



To zero-th order, the largest eigenvalue of L\ is given by the 
maximum of the function —r(v), whose expression is written 
in Eq. ©. The maximum of this function occurs at v = 0, 
and 



(o) _ _ 



= - r(0) " " n 



(36) 



The eigenfunction associated with this eigenvalue is indeed a 
delta function, centered in v = 0. In order to get the first 
order correction to the eigenvalue (f36t one has to project the 
zero-th order eigenfunction on the operator proportional to the 
small parameter: 

M « - J &viL x 5{vx) = ^= . (37) 

Finally one finds that for large A, the largest eigenvalue /z be- 
haves like: 

M (A) = -^(e- A -l) + 0(e- 2A ). (38) 

Hence, the probability of having Af collisions behaves, for 
values of Af small with respect to its average (Af) — ojt, as a 
Poisson distribution with a frequency equal tow/ 



P(Af) 



e ^2 / wt \ 



for TV < ut 



(39) 



For such a distribution, the large deviation function ir easily 
follows: 



7r(n) = n — nlog(nv2/iu) — Lo/y2 



(40) 



2. An approximate solution 

In order to get an approximate expression for the largest 
eigenvalue /i(A), we will suppose the associated eigenvector 
/(v, A) to be a Gaussian with a given temperature T(A). We 
expect that this approximation will provide accurate results 
for small values of A, given that for A = 0, /(v, 0) is exactly 
a Gaussian with a temperature To, but it is not a priori a sys- 
tematic approximation for larger values of A. Projecting the 
Boltzmann-like equation d28l ) onto the first two velocity mo- 
ments of the eigenfunction we want to compute, we are left 
with two closed equations for /z(A) and for T(A): 



/Lt(A) = u 



fi(X)dT(X) = v 2 



(41) 



where the v n are the collisional moments (the expression of 
the first ones is given in Appendix|E]i: 



(42) 



v n = J dvi / dv 2 / do-w"(vi 2 • <r)x 

/>r,A)^(vr)-/(vi,A)^(v 2 ) 

Solving simultaneously Eqs. ( |4TI ). one obtains: 



M(A) 



T2 y 



T(A) 




V2T 

vT+P 



(43a) 



(43b) 



One may notice that this result satisfies some of the previous 
requirements obtained from the asymptotic large A analysis. 
In particular /i(oo) = —uj/\f2, as in the previous subsection. 
Moreover, the fictitious temperature T(A) vanishes for infinite 
A, meaning that the eigenfunction associated with /j, does in- 
deed tend towards a delta function. Nevertheless, while these 
features concerning the zero-th order perturbation results of 
the previous section are fulfilled, the behavior of Eq. ( l43l is 
different at the next order. In fact, one can see from Eq. d43b 
that, for A — > oo: 



LO 



0(e 



(44) 



to be compared with Eq. (138b - This is a deficiency of the 
Gaussian approximation. 

In order to get an expression of w(n) one should compute 
the Legendre Transform of /i(A). This seems not feasible ana- 
lytically, but can be achieved numerically. The function n(n) 
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Jt(«)=min (ii,(X)+^, n) 
jt(h) uncorrelated (Eq. 23) 
Poisson(l) 
Poisson(l/V2) 



FIG. 3: (Color online) Plot of 7r(n) obtained in several fashions (a; = 
1). The solid line has been obtained from the numerical Legendre 
Transform of ^i(A) from Eq. <43t . The dashed line corresponds to 
the uncorrelated result obtained in Eq. d23 1 . while the dot-dashed line 
is the large deviation function of the Poisson distribution of average 
1. For completeness, we also show the large deviation function of 
the Poisson distribution with average l/y/2 (dotted line). 



is plotted in Fig. [3] together with the result of uncorrelated 
estimation of Eq. d23l l. as well as two Poisson distributions, 
of average to and u> / \2. One can see that the two estimations 
carried out are in very close agreement, apart from a slight 
difference in the behavior of the right tail. Moreover, the two 
new results presented here are clearly different from the Pois- 
son distribution with mean oj [denoted Poisson(l) since ui = 1 
in Fig. JQ], which is narrower, and underestimates extreme 
events, characterized by either very few or many collisions 
with respect to typical realizations. 



3. Sonine perturbation and cumulants 

We now further exploit the property that for A = the so- 
lution of the Boltzmann equation d28l l is exactly a Gaussian. 
This incites us, for A close to 0, to search for a solution /(v, A) 
as a small perturbation of a Gaussian distribution. One of the 
most useful expansions in kinetic theory is the Sonine poly- 
nomials expansion . In practice this expansion consists in 
looking for solutions expressed as a Gaussian times a series 
of Sonine polynomials, denoted S n (x): 



(2 \ 00 
' n=0 

The first Sonine polynomials are: 

S (x) = 1 , 



2TJ 



(45) 
(46a) 



Si(x) = ~ x +2 



S 2 (x) = -x 2 



2, d(d + 2) 



(46b) 



(46c) 



These polynomials have the property of being orthogonal with 
respect to a Gaussian measure in dimension d, and are there- 
fore related to Laguerre polynomials. From this feature, it 
follows that the coefficient of the series d45b do = 1 and that 
di = 0. Hence, the first nontrivial correction to the Gaussian 
approximation comes from the term proportional to 0,2 in the 
expansion d45l >. The procedure to get an estimate of the co- 
efficient a2 consists of solving a closed system of equations 
obtained projecting the equation d28l > onto the first velocity 
moments : 



fj,(X)m r , 



where: 



m„(A) 



dv/(v,A) , 



(47) 



(48) 



and v n (X) denotes the collisional moment of order n defined 
in Eq. d42b . Truncating the expansions d45b up to the second 
Sonine polynomial, one gets for the moments m n : 



m 



1 



TO2 = cff(A) , 



m 4 = d(d + 2)T 2 {\) (l + a 2 (A)) 



(49a) 



(49b) 



The expression of the first collisional moments in the Sonine 
approximation is given in the Appendix[E] Hence, taking the 
moment equation d48b for n = 0, 2 and 4 gives a closed sys- 
tem of equations for /i(A), T(X) and 0.2(A). This system can 
be solved perturbatively expanding its solutions in power se- 
ries around A = 0. The cumulants of JV obtained by means of 
the latter expansion turn out to be remarkably accurate. The 
expansion of fi(X) up to the 3 rd order in A gives access to the 
first three cumulants, which read: 



uit 
(AA 2 ) C 

cot 
(Af% 

tot 



1 , 
9 

64 



(50a) 



1 



4d + 3 / 
28d(64<i(320d + 729) 



(50b) 



35775) + 257391 



8192(4d+3) 3 



(50c) 



As can be noted from these values, the variance of the number 
of collisions already deviates by more than 12.5% from its 
Poisson value (equal to unity). 



IV. NUMERICAL RESULTS 



' For a general treatment of Sonine polynomials in kinetic theory, see e.g. 
|11] and references therein. 



In this section, we compare the theoretical results obtained 
in the previous sections against numerical simulations. In 
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FIG. 4: (Color online) The circles are the numerical measurements 
of At(A) for the Maxwell model. The solid line is the generating 
function of the Poisson distribution. The inset shows the same data, 
shifted vertically by one (in such a way that /i(A) + 1 is always pos- 
itive), in semi-logarithmic scale. 

addition to the two aforementioned techniques of Molecular 
Dynamics and Direct Simulations Monte Carlo (see Fig. [TJ, 
we have shaped a third numerical tool, constructing a Monte 
Carlo algorithm in order to directly solve the eigenvalue equa- 
tion (|28V We start by setting the stage for the latter method, 
before briefly commenting on the Molecular Dynamics sim- 
ulations used to measure the statistics of the number of colli- 
sions suffered by tagged particles. 

A. Monte Carlo approach 



In order to derive an algorithm for solving Eq. (l28l ). it is 
useful to rewrite this equation in the form: 

M (A)/(v, A) - e- A /[/>] + (1 - e- A )A[/>] , (51) 

where 

I[f\g] = J dv J d<r(v 12 • <x)[/r<?2* - fito] > (52) 

is the collision integral describing the elastic collision be- 
tween two particles having respectively a velocity pdf / and 
g. We are using the short-hand notation /** = f(v**). The 
functional is the loss term of the above collision inte- 

gral. It actually describes the statistics of hard spheres which 
annihilate after each collision: 

A[f\g] = -Jdvj d^(v 12 • &)frg 2 . (53) 

In the context of the nonlinear Boltzmann Equation, i.e. when 
4> = / in Eq. dBTt . and for positive values of A, the rhs of 
Eq. (|5"TT i exactly describes the time evolution of the veloc- 
ity pdf of the probabilistic ballistic annihilation process. This 
model, introduced in [3], consists in a system of N particles 




FIG. 5: (Color online) Numerical measurement of /i(A) for the hard- 
sphere model. The circles show the results of Monte Carlo simu- 
lations. The solid line is the theoretical prediction within the Gaus- 
sian approximation. The dashed line is the theoretical prediction (ob- 
tained solving numerically d47|48l > in the framework of the first So- 
nine correction). Finally the dotted line is the generating function of 
a Poisson distribution of parameter ui/\f2, which should asymptoti- 
cally dominate for large A (see Fig- [6j. The inset shows a zoom near 
A = 0, where the difference between Gaussian and Sonine orders 
can only be appreciated for the larger values of A displayed. 
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X 

FIG. 6: (Color online) Same data as in Fig. [5] divided by the gener- 
ating function of a Poisson distribution of parameter co/y/2. 



which move ballistically, and interact when at contact. The in- 
teraction may be an elastic collision, with probability e~ A , or 
an annihilation (with probability 1 — e~ A ). In the context of 
the linear Boltzmann equation, this process can be extended 
in the following way. Consider a set of N independent sys- 
tems, one of each is just a single particle, characterized by 
its velocity v%(t) (i — 1, . . . , N), assuming spatial homo- 
geneity. If each of these systems (particles) evolves in a hard 
sphere gas (in the thermodynamic limit) thermalized at tem- 
perature To, with a one point velocity pdf <j), and can both 
collide elastically (with probability e~ A ), or annihilate (with 
probability 1 — e~ A ), then, the reduced one point velocity pdf 
g(v, t) — (S^5(v — Vj(t))) will verify the Boltzmann Equa- 
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FIG. 7: (Color online) Temperature (defined as the variance) of the 
eigenfunction /(v,A) measured in Monte Carlo simulations. The 
solid and dashed lines are the predictions of the Gaussian and Sonine 
approximations respectively. 
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FIG. 8: (Color online) Coefficient a 2 (A-dependent velocity kurto- 
sis), measured in a Monte Carlo simulation. The dashed line corre- 
sponds to the numerical solution of the system l |49l >, which is correct 
only for small values of lo- 



tion 



dg{y,t) 
dt 



= L\g(v,t) 



(54) 



Since N(t) — J dvg(v, t) is the number of particles at time t, 
and since we know that for long times <?(v, t) ~ e M *, it is clear 
that for long times one has that N(t) ~ A^e^*. Note that the 
present interpretation implicitly assumes positive values of A. 
Moreover, the particles can only annihilate or collide; then the 
total number of particles can only decrease, and hence fi must 
be negative. The above observations provide a numerical tool 
for measuring /tt(A) as the decay rate of the total number of 
particles. The main difficulty, which we have successfully ad- 
dressed, with the above algorithm, is that the number of parti- 
cles constantly decreases, and there is no steady state but the 
trivial state N = 0. To circumvent this practical difficulty 
(that would lead to somewhat noisy statistics in the simula- 



tions), we have introduced an external source of particles (sys- 
tems), acting in such a way that the total number of particles 
is conserved. If every time that an annihilation takes place, a 
new particle is inserted as the clone of one of the N — 1 re- 
maining particles (chosen uniformly among this population), 
then the evolution equation of g(v, t) reads: 



<9t#(v, t) = £a5(v, t) + sg(v, t) 



(55) 



where s is a constant rate. At late times, if the above equation 
has a steady state, the largest eigenvalue of the operator L \ + s 
vanishes, and hence one has that /i(A) + s = 0. Finally one 
can measure /i(A) simply as (minus) the steady state average 
of the number of particles injected by the external source. 

So far, we showed how to construct a Markov chain in order 
to simulate the eigenvalue equation d28b for positive values of 
A. For negative values of A the procedure is almost identical. 
Introducing a new time scale t — (2e _A — l)t, Eq. d54l i can 
be rewritten as: 



2e~ A - 1' 



1 



2e- 



— AW] 



(56) 



All the previous physical interpretations and remarks still 
hold, except that now instead of having annihilation (with 
probability (e~ A - l)/(2e~ A - 1)), one has duplication (or 
cloning). Hence /i(A) will be positive, and one can add an ex- 
ternal source in order to remove particles when new particles 
are created. Summarizing, the algorithm proceeds as follows: 

(0) The velocity of the N particles are stored in a JV x d 
matrix. A scalar s is set to 0. A (small) time step dt is 
chosen. 

(1) A particle is chosen randomly in the population with 
uniform probability. Its velocity is denoted Vi . 

(ii) An interaction is accepted with a probability oc v 12 • 
<rO(vi2 • cr) x dt, where <r is a random direction in d 
dimensions, and v 2 a d-dimensional zero-mean Gaus- 
sian random variable of variance To. In practice, dt has 
to be chosen in such a way that |vi2 ■ cr\dt is always 
smaller than (or equal to) one. 

(iii) When the interaction is accepted, if A > (resp. 
A < 0), the particle will have a post-collisional veloc- 

with probability e~ A ^resp. 9p I A _ 1 ) ) or will be 



ltyvj „x^^„„^^ V w 2e _ 
removed (resp. duplicated) otherwise. 

(iv) If the particle has been removed (resp. duplicated) in 
step (iii), one of the N — 1 remaining particles is cho- 
sen randomly and uniformly, and is duplicated (resp. 
removed), s is increased by 1 . 

(v) Time is increased by an amount dt (resp. dt). 

(vi) Back to (i). 

This algorithm bears some similarities with the approach pro- 
posed in IU2I1 . also intended to directly measure large devi- 
ation functions. Nonetheless, the version proposed here is 
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more inspired by some variants of the DSMC algorithm for 
systems which do not conserve the total number of parti- 
cles QHEl. In order to check the reliability of the al- 
gorithm, we have first performed simulations in the case of 
Maxwell molecules, where the number of collisions is exactly 
distributed following the Poisson distribution (cf. the Intro- 
duction). The measurements of /i(A) for this particular model 
are shown in Fig. HI together with the generating function of 
the Poisson distribution. The agreement is excellent. In the 
simulations the temperature scale is set by the temperature of 
the heat bath T , which we set to unity. The time scale is set 
by the mean free time, which we also set to unity. In the sim- 
ulation data, the mean collision frequency is therefore u = 1. 

In the case of the hard-sphere model, a measurement of the 
largest eigenvalue /i(A) is shown in Fig. [5] One can see that 
the numerical results are in very good agreement both with 
the Gaussian and the Sonine approximations. Figure [6] shows 
the large A behavior of ^t(A) divided by the prediction ( |38T >. 
In this case, a good agreement between the numerical data 
and the theoretical predictions is found. In the framework of 
the above described Monte Carlo algorithm, it is also possible 
to measure the stationary velocity pdf g(v, t — oo, A) which 
is equal to the eigenfunction /(v, A) associated with /i(A). 
Hence one can also compare the analytical predictions for the 
temperature T(A) and for the Sonine coefficient 02(A) with 
the Monte Carlo results, see Figs. [7]and[8] In Fig. [7] one can 
see the temperature T(A). Here the Gaussian approximation is 
already able to capture the general behavior of this "effective" 
temperature, and for small values of A, the Sonine corrections 
compare very well with the simulation data. In Fig. [8] one 
can see the coefficient 02 as a function of A. As expected, 
for small values of A the expansion carried out in the previous 
section correctly describes the result of the simulations. 
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FIG. 9: (Color online) Large deviation function 7r(n) of the number 
of collisions M suffered by each single particle. Symbols show the 
results of Molecular Dynamics simulations. The solid line is the re- 
sult in the framework of the Gaussian approximation, and the dashed 
line is the large deviation function of the Poisson distribution, of av- 
erage 1: 7r(n) = n — nlogn — 1. 



TABLE I: Cumulants for the number of collisions M from MD sim- 
ulations, and comparison with the Gaussian and Sonine approxima- 
tions (time is measured in units of the mean-free time). 





(N) c /t 


(N 2 ) c /t 




t = 10 


1. 


1.1228 


1.1282 


t = 50 


1. 


1.1354 


1.1045 


Poisson 


1 


1 


1 


Gaussian 


1 


1.125 


1.1289 


Sonine 


1 


1.1377 


1.1073 



B. Molecular Dynamics 

The most direct way to measure the large deviation func- 
tion of the number of collisions JV is of course to count it in 
a (numerical) experiment, and then construct its probability 
distribution function. To this end, we have performed Molec- 
ular Dynamics (MD) simulations of a two dimensional hard- 
disk gas of A = 10 3 particles of diameter a = 1 at density 
p = N/V = 1CT 3 <T~ 2 and 10~ 2 cr~ 2 . The particles evolve 
in a square box with periodic boundary conditions, and the 
time is measured in mean free time units, in such a way that 
(JV) = t. We measured the statistics of the number of colli- 
sions JV suffered by each particle. The value of the first cu- 
mulants is reported in Table [I] together with the Poisson pre- 
diction, as well as the results from the Gaussian and Sonine 
approximations. It seems that, when time increases, the cu- 
mulants converge towards a value which is close to the Sonine 
predictions. However, it must be noted that even if increasing 
time makes finite time corrections smaller, the statistics be- 
come poorer and poorer. The agreement with the results from 
the Sonine approximation is very good, while the Gaussian 
order already provides a reliable estimation. This is further 
shown in Figure |9j where two measurements of the large de- 



viation function of the number of collisions for two different 
times, t = 10 and t = 50 (with a = 1 and To = 1), are com- 
pared to the numerical inverse Legendre Transform of /i(A) at 
Gaussian order (given by Eq. d43l)). Here again, the Poisson 
prediction is distinctly off. 



V. CONCLUSION 

We have shown that the collisional statistics of the hard 
sphere gas exhibits clear deviations from the Poisson distri- 
bution. These deviations have been consistently quantified 
both analytically and numerically. In the analytical treatment, 
the cumulant generating function /i(A) plays a pivotal role. 
The eigenvalue equation defining this object can be fruitfully 
interpreted in terms of population dynamics with annihila- 
tion and cloning events, which eventually leads to an effi- 
cient algorithm allowing to compute the various quantities in- 
volved in the theoretical analysis. The corresponding numeri- 
cal method, of Monte Carlo type, should not be confused with 
the more conventional Direct Simulation Monte Carlo tech- 
nique, which we also implemented, and that is intended to 
solve a different kinetic equation (the Boltzmann equation). 
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Finally, a third numerical method (Molecular Dynamics) was 
employed. These three routes provide complementary and 
valuable results, that strongly support our predictions. 

Interestingly, the present formalism can be extended to the 
study of out of equilibrium systems, such as granular gases, 
where the question of the collisional statistics has been the 
focus of recent interest G3 H d El HI. In particu- 
lar, the strong effect of dissipation on the distribution of free 
flight times reported in lll9l l20ll calls for further investiga- 
tions. In addition, if one considers a gas of inelastic smooth 
hard spheres, kept in a steady state by a velocity independent 
force (as e.g. a vibrating wall of the container, or a stochastic 
force acting independently on each particle), then the phase 
space volume of the system has been reduced, after a time t, 
by a factor (1 — a)^^\ where a is the coefficient of normal 
restitution 11611 . Hence one sees that in this non-equilibrium 
system, the number of collisions can be exactly identified, up 
to a constant prefactor, with the integrated phase space con- 
traction rate. This quantity has already been the subject of 
many works in non-equilibrium statistical mechanics, and is 
often intimately related with the irreversible entropy produc- 
tion (see e.g. fl22ll and references therein). 
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and ([Tol l, yields 

dt t / dv e- r W 0(v) = -L, (A3) 

o J u (r) 

which is the required relation. 

Second, expression ( IA1I ) explicitly slightly differs from the 
correct free flight time distribution in two limiting cases. 

• At small values of t, Pppx behaves as: 

P¥ FT it) ~ (r(v)) ~ (r(vf)t + 0{t 2 ) , (A4) 
while the true distribution behaves as 



Pfft (t) 



to 



(A5) 



so that even the value in t — is different, although this 
difference is numerically small (see also Fig. ITOb . For 
instance, in d = 2 one has: 



(r(v) 2 ) 



and 



(r(v) 3 ) 
to(r(v) 2 



1.06354 



1.13962 



(A6) 



(A7) 
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APPENDIX A: COMMENT ON THE FREE FLIGHT TIME 
DISTRIBUTION 

We shall give here some additional arguments on the in- 
correctness of Eq. (O, when used with a Gaussian weight </>, 
instead of the weight <fi co u defined in Eq.lfTOl). We shall refer 
to this (erroneous) distribution as Pppx- 



Pfft (0 



dv</>(v)P(i|v) 



(Al) 



A first argument bears on the inconsistency between the above 
relation and the definition of the collision frequency co = 
(r(v)}, where the brackets denote an average over a Gaussian 
weight. Indeed, the average time r between two subsequent 
collisions (mean free time) of a given particle is equal to the 
inverse of the collision frequency: r = l/(r(v)). Besides, 
from Eq. ( lAlb the mean free time r is obtained as: 



At tP, 



w 



FFT 



(*) = 



1 



1 



r{v)l (r(v)) 



(A2) 



Hence we see that expression (IAU is in contradiction with the 
definition of the collision frequency. On the other hand, the 
counterpart of Eq. ( IA2b with the distribution provided by (0 



• The large time behavior of the probability also slightly 
differs from Eq. ( ITTb : 



Pfft ~ ex P 



Lut\ d-i ( 2 y/2wt\ 



-d/2 



(A8) 

In the above expression, the leading exponential term 
remains the same as in Eq. ( fTTT i. since it is only deter- 
mined by the minimum of the function r(v). Nonethe- 
less, the subleading prefactor is slightly different, as it 
can be appreciated in Fig. [TOl 



APPENDIX B: FREE FLIGHTS TIME DISTRD3UTION FOR 
VERY HARD PARTICLES 

For the sake of completeness, we report in this appendix 
the result for the free flight time distribution arising in the 
very hard particle (VHP) model. This framework allows for 
an explicit analytic computation of Pfft^)- The model, in- 
troduced by Ernst and Hendriks 112 311 . consists in a choice of 
the collision rate proportional to the total energy of the sys- 
tem. In the case of a single tagged particle, the corresponding 
linear Boltzmann equation reads: 

|/(vi,t) = \ |dv 2 Ja& (vi ^ )2 [f?tfr - hh] . 

(Bl) 

The velocity dependent collision rate is thus quadratic in v: 



rvHp(v) _ 1 v 2 
(^vhp 2 2effb 



(B2) 
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FIG. 10: (Color online) Free flights time distribution of a hard disc 
gas. The circles corresponds to the results of Molecular Dynamics 
(MD) simulations at density p — 0.01<j~ 2 . The full line is the nu- 
merical integration of Eq. (J9}. The dashed line is the numerical inte- 
gration of Ppp T , and the dotted line its large time behavior described 
by Eq. ( IA8b . The inset shows the same data in semi-logarithmic 
scale, where it is possible to note that at large times the prediction 
of Ppft is always above the true distribution. Time is measured in 
units of the mean free time (u) — 1). 



where 



WVHP 



£T(d/2) 



(B3) 



is the collision frequency for VHP. Due to the simpler form 
of the collision rate, the free flight time distribution can be 
expressed analytically: 



^VHP 
FFT 



(t) = -d^e-^x 
y ' 4 

{d + t)~i~ 2 (Ad 2 + (4t + 2)d + t 2 ) , (B4) 



where the time is here in units of the mean collision rate 
luvhp- Note that here the large time behavior is ~ e -^*/ 2 ; 
which is even slower than the Hard-Sphere gas behaviour 
(which is ~ e -"*/v2) w ith respect to the Poissonian case 
(~ e~ wt ). A comparison between these three distributions 
is shown in Fig. QT| The observed large time behaviors sug- 
gest that Maxwell and VHP models provide a upper and lower 
bounds for the hard-sphere model, a phenomenon reminiscent 
to that observed in the long time behaviour of dynamics that 
do not conserve the density lfl5ll . 



APPENDIX C: SADDLE-POINT APPROXIMATION FOR 
THE COLLISION RATE 

In this appendix we show how to recover Eq. (TT~2T > with the 
saddle point method. We first note that r(v) can be expressed 
as the difference of two different integrals: 



r{v) =n d _ 1 {I 1 (v)-I 2 {v)) 



(CI) 



- Very Hard Particles 

- Hard Spheres 

■■ Maxwell Molecules 




FIG. 1 1 : (Color online) Free flights time distributions for three dif- 
ferent interaction kernels, in two dimensions. The (red) dashed line 
shows the result for VHP (Eq. IB4b . the (blue) solid line and the 
(green) dotted line show respectively the results for Hard Spheres 
and Maxwell Molecules. Time is measured in units of the mean free 
time {uj — 1). 



where 



h{v) 



/ + 00 />TT 
dv 2 a / d0sin0 d - 2 (ucos0)x 
-oo JO 

6(VC0S6) - V 2<J )(j)(v 2(7 ) , (C2) 



h(v) = 



dv 2a - I ddsm9 d 2 v 2<J x 
'o 



0(UCOS6» - V 2a )4>{ v 2a) ■ (C3) 



For the first integral it is more convenient to perform first the 
integration over the angle 6, which leads to: 



h(v) 



d-1 



+ °° dv 2a (\~(^L) 2 \ 2 0( U2ff ) 



(C4) 

Then, defining the rescaled variables v = v/^/2Tod and 
v 2rT = v 2a I V2Tb<i one finds, to leading order in d: 



h (v) ~ v 



+oo 



exp 



dv 2a x 




(C5) 



When d is very large, the above integral is dominated by the 
maximum of the function inside the exponential, which turns 
out to be located in v 2a = 0. One can then perform a series 
expansion around ii 2a = up to the second order. This results 
in a Gaussian integral which is easily integrated, and yields: 



h(v) 



I To 2v 



-2 



d vl + 2v 2 



(C6) 
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FIG. 12: (Color online) Free path length distribution for Maxwell 
molecules in two dimensions. The circles are the result of a DSMC 
simulation, while the solid line is the numerical integration of Eq. 
l |D3l >, The dashed line is the large length approximation of Eq. l lD4t . 
All data are reported in units of mean free path. 



As for the second integral involved in the expression of r(v) 
the simplest is first to perform the integration over v^; hence 
one obtains, to leading order in d: 

I 2 (v) = -\ ^ [ exp [eZ(-u 2 cos 2 + lnsin0)] . 

V 27T Jo 

(C7) 

The only maximum of the function in the exponential between 
and 7T is in 9 = ir/2. Then, expanding as usual this function 
around the maximum up to second order, and extending the 
range of integration from — oo to +oo, one finds: 



1 



(C8) 



and the free path length distribution is simply the average of 
the above probability over a Gaussian weight (for Maxwell 
molecules the on collision distribution <\> co u is still a Gaus- 
sian): 



Pfpl(x) = / d^(v)pV pL (x\v) . (D3) 



This last expression can be expressed analytically in terms of 
Meijer G-functions [7], but here we shall focus only on the 
large length behavior, for which simpler expressions are avail- 
able. In particular, the saddle point approximation gives: 



pM 



FPL 




(D4) 



which leads to a stretched exponential behavior at large x (s 
cxp-a; 2 / 3 )(cfFig.[l2|. 



APPENDIX E: COLLISIONAL MOMENTS 

Here we provide the expressions of the first collisional mo- 
ments, defined by Eq. ( t42t . for / a Gaussian, and a Gaussian 
multiplied by a Sonine Polynomial. We denote by 



dvi<i A 



(27rT) d / 2 



(El) 



Finally, noting that T(d/2)/T((d - l)/2) ~ v/f.Eq. O 
follows. 



{2irT) d / 2 



2T J 



(E2) 



APPENDIX D: FREE PATH LENGTH DISTRIBUTION FOR 
MAXWELL MOLECULES 

In this appendix, we investigate the free path length distri- 
bution for Maxwell molecules. As already mentioned, in this 
case the collision rate is a constant, cj, independent of the ve- 
locity of the particle. It follows then that the free flight time 
probability is exponential: 



pM 



FFT 



pM 



FFT 



(t) 



(Dl) 



where L\ is the operator appearing in d28l i. Hence in the 
Gaussian approximation one has v n = , while in the So- 
nine approximation v n = + a-i • The expressions of 

(i) 

the first Vn are: 

„(0) _ I ( e ~ 



i VT + n 



'2-K 



(E3) 



Besides, the free path length distribution for a given velocity 
v reads: 

P F 4 PL (x\v) = ~exp( c ^x) , (D2) 



CO = r 2 + 2Tr + 2r 2 ~e A r (3T + 2T ) 

' 2 ~ e A V2~7Ty/T + T a 



(E4) 
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vf ] = j-y/Tr (- (e A T 2 (8T 2 + 12 TT + 3 T 2 )) + T 2 (3 T 2 + 12 TT + 8 T 2 )) - 

4T(T + To) / 3(-i + e A ) y¥r(r + r ) | (e A r (2T + r ) -r (r + 2T )) ' 



\/2e A 7r(T + T ) i J 



(E5) 



v o — 



(-l + e A ) T 2 
8e A V2^{T + T ) 



(1) _ T 2 (-3 (l + e A + 2 (-l + e A )) T 2 - (-6 + 22e A ) TT -2 (-3 + 8e A ) T 2 ) 

8e A V27(T + T )5 



(E6) 



(E7) 



^ = |t 2 (-45 (-l + 5e A ) T 4 + 6 (28 - 132 e A ) T 3 T - (-228 + 1000e A ) T 2 T 2 + 4 (36 - 128 e A ) TT 3 + 

+8 (3-8e A ) T 4 ) J j |8e A V2^(T + T ) 5 | (E8) 
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